This document provides the code for the analyses reported in:

Cosme et al. (Preprint) Testing the adolescent social reorientation model using hierarchical growth curve modeling with parcellated fMRI data

This script reproduces the analyses reported in supplementary material.

load packages

library(tidyverse)
library(knitr)
library(lme4)
library(lmerTest)

define color palettes and labels

pal_self_other = c("#FFA90A", "#247BA0")
pal_social_academic = c("#63647E", "#F25F5C")
pal_wave = c("#693668", "#A74482", "#F84AA7")
pal_label = c("#47A8BD", "#DBC057", "#FF3366")
pal_gender = c("#70c1b3","#247BA0")

parcel_labeller = labeller(label = c('social' = 'social parcels', 'other' = 'control parcels', 'self' = 'self parcels'),
                           domain = c('social' = 'social domain', 'academic' = 'academic domain'),
                           wave = c("t1" = "wave 1", "t2" = "wave 2", "t3" = "wave 3"))

dcbw = theme_classic() +
  theme(text = element_text(size = 14, family = "Futura Medium", color = "black"),
        panel.background = element_blank(),
        plot.background = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.background = element_rect(fill = NA, color = NA),
        axis.line = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        panel.grid.minor = element_blank())

load MRI data

  • exclude participants with >20% volumes with motion artifacts or low quality FX data
  • standardize betas within parcel (across participants and time)
  • recode standardized betas +/- 3 SDs from the mean as NA (0.8% of all observations)
# define parcels
self_parcels = c(5, 17, 28, 47, 52, 66, 116, 147, 152, 156, 184, 198, 207, 225, 249, 292, 309, 354, 380)
social_parcels = c(18, 23, 29, 49, 54, 59, 62, 63, 67, 76, 111, 114, 139, 143, 146, 150, 178, 179, 189, 203, 212, 224, 229, 236, 238, 239, 245, 250, 259, 266, 271, 301, 305, 310, 322, 324, 328, 331, 333, 342, 343, 350, 374, 391)

# mri exclusions
mri_exclusions = c('s002_t1', 's004_t1', 's008_t1', 's011_t1', 's017_t1', 
                   's026_t1', 's033_t2', 's034_t1', 's041_t1', 's044_t1', 
                   's047_t1', 's051_t1', 's054_t1', 's057_t1', 's059_t1', 
                   's061_t1', 's063_t1', 's070_t2', 's074_t1', 's074_t2', 
                   's078_t1', 's084_t1', 's090_t2', 's090_t3', 's094_t1', 
                   's094_t2', 's096_t1') 

# mri exclusions
parcellations = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other'))) %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

# exclude parameter estimates 3 SD from the mean
parcellations_ex = parcellations %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

load cleaned SPPC and task data

# demographics
race_ethnicity = readxl::read_xlsx("../data/Ethnicity_w1.xlsx") %>%
  extract(SID, "subjectID", "L([0-9]{3})-.*") %>%
  mutate(subjectID = sprintf("s%s", subjectID),
         hispanic_latinx = ifelse(Ethnicity == 1, "yes", "no"),
         race = recode(Race, `1` = "White", `2` = "Multiracial", `3` = "Pacific Islander",
                       `4` = "Native American", `5` = "Black or African American")) %>%
  select(subjectID, hispanic_latinx, race)

demo = read.csv("../data/SFIC_age.pds.gender.csv") %>%
  rename("subjectID" = SID,
         "wave" = wavenum,
         "gender" = Gender) %>%
  mutate(subjectID = sprintf("s%03d", subjectID),
         wave = paste0("t", wave),
         age_c = age - 13,
         age_c2 = age_c^2,
         pdss_c = pdss - 3,
         pdss_c2 = pdss_c^2) %>%
  full_join(., race_ethnicity)

# self-report
sppc = readRDS("../data/sppc.RDS")

# task
task = readRDS("../data/task.RDS") %>%
  ungroup() %>%
  mutate(percent_stdgm = (task_percent - mean(task_percent, na.rm = TRUE)) / 10, #center at grand mean, in units of 10%
         percent_std50 = (task_percent - 50) / 10,
         percent_std = (task_percent - 50) / 10) %>% #center at 50, in units of 10%
  group_by(target, domain) %>%
  mutate(percent_stdc = (task_percent - mean(task_percent, na.rm = TRUE)) / 10) #center at condition mean, in units of 10%

# merge data
merged = parcellations_ex %>%
  full_join(., task, by = c("subjectID", "wave", "domain", "target")) %>%
  full_join(., sppc, by = c("subjectID", "wave", "domain")) %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) #no MRI, task, or self-report data was collected

# subset data for modeling
neuro_model_data_pdss = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex_pdss = neuro_model_data_pdss  %>%
  na.omit()

comp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, competence) %>%
  unique() %>%
  na.omit()

imp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, importance) %>%
  unique() %>%
  na.omit()

status_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, target, domain, task_percent, contains("percent")) %>%
  unique() %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy_pdss = neuro_model_data_ex_pdss %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

sample descriptives

full sample

merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave) %>%
  summarize(n = n(),
            meanAge = mean(age, na.rm = TRUE),
            sdAge = sd(age, na.rm = TRUE)) %>%
  kable(format = 'pandoc', digits = 2)
wave n meanAge sdAge
t1 90 10.07 0.31
t2 57 13.04 0.31
t3 44 16.46 0.58
merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave, gender) %>%
  summarize(n = n()) %>%
  spread(gender, n) %>%
  kable(format = 'pandoc', digits = 2)
wave F M
t1 45 45
t2 30 27
t3 25 19
merged %>% 
  select(subjectID, hispanic_latinx, race) %>%
  unique() %>%
  group_by(hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)
race no yes
Black or African American 3.3 NA
Multiracial 11.1 5.6
Native American 4.4 1.1
Pacific Islander 4.4 1.1
White 43.3 25.6

broken down by inclusion

age_table = merged %>% 
  select(subjectID, inclusion, wave, age, pdss, gender) %>%
  unique() %>%
  group_by(inclusion, wave) %>%
  summarize(N = n(),
            Mage = mean(age, na.rm = TRUE),
            SDage = sd(age, na.rm = TRUE),
            Mpdss = mean(pdss, na.rm = TRUE),
            SDpdss = sd(pdss, na.rm = TRUE))

gender_table = merged %>% 
  select(subjectID, inclusion, wave, age, pds, gender) %>%
  unique() %>%
  group_by(inclusion, wave, gender) %>%
  summarize(N = n()) %>%
  spread(gender, N) %>%
  rename("Female" = F,
         "Male" = M)

merged %>% 
  select(subjectID, inclusion, hispanic_latinx, race) %>%
  unique() %>%
  group_by(inclusion, hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)
inclusion race no yes
completed MRI Black or African American 0.8 NA
completed MRI Multiracial 7.6 2.5
completed MRI Native American 3.4 0.8
completed MRI Pacific Islander 2.5 0.8
completed MRI White 27.7 16.0
didn’t complete MRI Black or African American 2.5 NA
didn’t complete MRI Multiracial 1.7 NA
didn’t complete MRI Native American NA 0.8
didn’t complete MRI Pacific Islander 0.8 NA
didn’t complete MRI White 8.4 3.4
excluded from MRI Multiracial 0.8 2.5
excluded from MRI Native American 0.8 NA
excluded from MRI Pacific Islander 0.8 0.8
excluded from MRI White 11.8 2.5
age_table %>%
  left_join(gender_table) %>%
  extract(wave, "Wave", "t([1-3]{1})") %>%
  arrange(Wave) %>%
  select(Wave, inclusion, N, Mage, SDage, Mpdss, SDpdss, Female, Male) %>%
  rename("MRI status" = inclusion) %>%
  kable(format = 'pandoc', digits = 2)
Wave MRI status N Mage SDage Mpdss SDpdss Female Male
1 completed MRI 57 10.08 0.32 2.16 0.74 33 24
1 didn’t complete MRI 12 10.00 0.25 2.00 0.77 4 8
1 excluded from MRI 21 10.07 0.33 1.83 0.58 8 13
2 completed MRI 44 13.06 0.33 3.57 1.05 25 19
2 didn’t complete MRI 8 13.07 0.20 3.69 0.65 4 4
2 excluded from MRI 5 12.79 0.28 2.62 0.85 1 4
3 completed MRI 34 16.33 0.46 4.75 0.39 19 15
3 didn’t complete MRI 9 17.06 0.60 4.78 0.67 6 3
3 excluded from MRI 1 15.72 NA 3.50 NA NA 1

age and puberty correlations

cor_fun = function(data) pmap(var.names, ~ cor.test(data[[.x]], data[[.y]])) %>% 
  map_df(broom::tidy) %>% 
  cbind(var.names, .)

age_pds = merged %>%
  select(subjectID, wave, gender, age, pdss) %>%
  gather(dev, value, age, pdss) %>%
  unite(wave_dev, dev, wave, sep = " ", remove = TRUE) %>%
  rownames_to_column() %>%
  spread(wave_dev, value) %>%
  group_by(subjectID) %>%
  fill(everything(), .direction = "up") %>%
  fill(everything(), .direction = "down") %>%
  select(-rowname) %>%
  unique() %>%
  ungroup()

var.names = tidystringdist::tidy_comb_all(names(select(age_pds, -c(subjectID, gender)))) %>%
  filter(!(grepl("age", V1) & grepl("age", V2)))

age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V2, r) %>%
  arrange(Sex) %>%
  filter(Sex == "F") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")
Sex V1 PDS-S 1 PDS-S 2 PDS-S 3
F Age 1 .50*** .20 .01
F Age 2 -.12 .40* .30
F Age 3 -.11 .02 .24
F PDS-S 1 –- .42* -.08
F PDS-S 2 –- –- .25
age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V1, r) %>%
  arrange(Sex) %>%
  filter(Sex == "M") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")
Sex V2 Age 1 Age 2 Age 3 PDS-S 1 PDS-S 2
M PDS-S 1 -.32* -.37 -.26 –- –-
M PDS-S 2 .22 .41* .50* .16 –-
M PDS-S 3 -.24 .37 .41 .30 .62**

plot

full sample

merged %>% 
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(19, 17, 21)) +
    labs(x = "\npubertal stage (PSD score)", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

MRI only

merged %>% 
  filter(!inclusion == "didn't complete MRI") %>%
  mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(21, 19)) +
    labs(x = "\npubertal stage (PSD score)", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

neural puberty models

visualize raw

Visualize the developmental trajectories using the raw data

hypothesis 0

domain_parc_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_raw = cowplot::plot_grid(domain_parc_plot_raw, target_parc_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, group = domain, color = domain)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_social_academic[2], fill = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h1_raw = cowplot::plot_grid(domain_plot_raw, soc_acad_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, group = target, color = target)) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, beta_std_avf) %>%
  unique() %>%
  spread(target, beta_std_avf) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_self_other[2], fill = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h2_raw = cowplot::plot_grid(target_plot_raw, self_other_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, beta_std, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = beta_std, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_y_continuous(breaks = c(-.4, .2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', linetype = '', fill = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_int_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, beta_std_avg) %>%
  unique() %>%
  spread(target, beta_std_avg) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_raw = cowplot::plot_grid(int_plot_raw, soc_acad_int_plot_raw, self_other_int_plot_raw,
                       labels = c('A', 'B', 'C'), ncol = 3))

run models

Estimate full models

domain x puberty

model_domain_equation = formula(beta_std ~ 1 + pdss_c*domain*label +
                                  pdss_c2*domain*label +
                                  (1 + pdss_c*domain || subjectID) +
                                  (1 + pdss_c*domain + pdss_c2*domain || parcellation))

model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy_pdss)

model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 10*(model_domain_numFx + model_domain_numRx + 1)^2

if (file.exists("../data/model_domain_pdss.RDS")) {
  model_domain = readRDS("../data/model_domain_pdss.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = neuro_model_data_dummy_pdss, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain_pdss.RDS")
}

domain x target x puberty

model_target_equation = formula(beta_std ~ 1 + pdss_c*target*domain*label +
                                  pdss_c2*target*domain*label +
                                  (1 + pdss_c*target*domain || subjectID) +
                                  (1 + pdss_c*target*domain + pdss_c2*target*domain || parcellation))

model_target_formula = lFormula(model_target_equation, data = neuro_model_data_dummy_pdss)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

if (file.exists("../data/model_target_pdss.RDS")) {
  model_target = readRDS("../data/model_target_pdss.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy_pdss, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target_pdss.RDS")
}

compare

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
model_domain 29 438603.3 – – –
model_target 57 438378.3 280.92 28 < .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (puberty 3, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("pdss_c", "Puberty", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (puberty 3, label (control)) -0.008 [-0.060, 0.045] 0.027 -0.280 398.452 .779
Puberty 0.006 [-0.005, 0.017] 0.005 1.122 114.748 .264
Target -0.012 [-0.030, 0.006] 0.009 -1.343 163.818 .181
Domain 0.004 [-0.021, 0.030] 0.013 0.342 126.671 .733
Label (self) -0.010 [-0.214, 0.194] 0.104 -0.097 356.861 .923
Label (social) 0.209 [0.069, 0.348] 0.071 2.936 356.762 .004
Puberty2 -0.000 [-0.006, 0.006] 0.003 -0.013 699.220 .990
Puberty x Target 0.004 [-0.006, 0.014] 0.005 0.833 87.273 .407
Puberty x Domain 0.004 [-0.013, 0.021] 0.009 0.478 65.176 .634
Target x Domain -0.009 [-0.046, 0.027] 0.019 -0.506 146.016 .614
Puberty x Label (self) 0.004 [-0.021, 0.029] 0.013 0.308 378.167 .759
Puberty x Label (social) 0.015 [-0.002, 0.033] 0.009 1.705 377.575 .089
Target x Label (self) 0.214 [0.157, 0.271] 0.029 7.340 2050.846 < .001
Target x Label (social) -0.019 [-0.058, 0.020] 0.020 -0.949 2021.793 .343
Domain x Label (self) 0.126 [0.058, 0.194] 0.035 3.612 692.867 < .001
Domain x Label (social) 0.144 [0.097, 0.191] 0.024 6.060 685.991 < .001
Target x Puberty2 -0.002 [-0.010, 0.006] 0.004 -0.508 360.135 .612
Domain x Puberty2 -0.012 [-0.022, -0.002] 0.005 -2.255 558.945 .025
Label (self) x Puberty2 -0.012 [-0.029, 0.004] 0.008 -1.480 383.964 .140
Label (social) x Puberty2 -0.001 [-0.012, 0.010] 0.006 -0.174 380.939 .862
Puberty x Target x Domain -0.019 [-0.041, 0.004] 0.012 -1.610 76.241 .112
Puberty x Target x Label (self) 0.027 [-0.002, 0.055] 0.015 1.816 176512.793 .069
Puberty x Target x Label (social) 0.008 [-0.011, 0.028] 0.010 0.839 176505.823 .402
Puberty x Domain x Label (self) 0.021 [-0.007, 0.050] 0.015 1.470 176512.731 .142
Puberty x Domain x Label (social) 0.007 [-0.013, 0.027] 0.010 0.699 176504.705 .484
Target x Domain x Label (self) -0.128 [-0.239, -0.018] 0.056 -2.280 176522.246 .023
Target x Domain x Label (social) 0.055 [-0.020, 0.130] 0.038 1.432 176510.597 .152
Target x Domain x Puberty2 0.005 [-0.012, 0.022] 0.009 0.600 381.300 .549
Target x Label (self) x Puberty2 -0.022 [-0.047, 0.004] 0.013 -1.647 176520.765 .100
Target x Label (social) x Puberty2 0.002 [-0.016, 0.019] 0.009 0.186 176510.160 .852
Domain x Label (self) x Puberty2 -0.003 [-0.029, 0.023] 0.013 -0.204 1002.405 .838
Domain x Label (social) x Puberty2 0.005 [-0.013, 0.022] 0.009 0.512 989.901 .608
Puberty x Target x Domain x Label (self) -0.041 [-0.099, 0.016] 0.029 -1.414 176512.990 .157
Puberty x Target x Domain x Label (social) -0.025 [-0.064, 0.014] 0.020 -1.263 176504.330 .207
Target x Domain x Label (self) x Puberty2 0.101 [0.050, 0.152] 0.026 3.849 176520.602 < .001
Target x Domain x Label (social) x Puberty2 0.020 [-0.015, 0.055] 0.018 1.123 176506.142 .261
# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)
##  Groups          Name                  Variance  Std.Dev.
##  parcellation    target:domain:pdss_c2 0.0000000 0.00000 
##  parcellation.1  pdss_c:target:domain  0.0000000 0.00000 
##  parcellation.2  domain:pdss_c2        0.0000992 0.00996 
##  parcellation.3  target:pdss_c2        0.0000000 0.00000 
##  parcellation.4  target:domain         0.0000000 0.00000 
##  parcellation.5  pdss_c:domain         0.0000000 0.00000 
##  parcellation.6  pdss_c:target         0.0000000 0.00000 
##  parcellation.7  pdss_c2               0.0004905 0.02215 
##  parcellation.8  domain                0.0075317 0.08679 
##  parcellation.9  target                0.0010029 0.03167 
##  parcellation.10 pdss_c                0.0020591 0.04538 
##  parcellation.11 (Intercept)           0.1899332 0.43581 
##  subjectID       pdss_c:target:domain  0.0029573 0.05438 
##  subjectID.1     target:domain         0.0054222 0.07364 
##  subjectID.2     pdss_c:domain         0.0027859 0.05278 
##  subjectID.3     pdss_c:target         0.0003843 0.01960 
##  subjectID.4     domain                0.0039597 0.06293 
##  subjectID.5     target                0.0010293 0.03208 
##  subjectID.6     pdss_c                0.0006930 0.02633 
##  subjectID.7     (Intercept)           0.0024120 0.04911 
##  Residual                              0.6669937 0.81670

visualize fitted models

Visualize the developmental trajectory using the fitted values from the domain x target x age model

reForm = as.formula("~(1 + pdss_c*target*domain + pdss_c2*target*domain || parcellation)")

neuro_plot_data_puberty = with(neuro_model_data_dummy_pdss, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                pdss = unique(pdss),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         pdss_c = pdss - 3, 
         pdss_c2 = pdss_c^2, 
         subjectID = NA)

neuro_plot_data_puberty$expected = predict(model_target, newdata = neuro_plot_data_puberty, re.form = reForm)
neuro_plot_data_puberty$expected_mean = predict(model_target, newdata = neuro_plot_data_puberty, re.form = NA)

neuro_plot_data_puberty = neuro_plot_data_puberty %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))

hypothesis 0

domain_parc_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = domain)) +
  #geom_point(aes(group = target)) + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  scale_color_manual(values = pal_social_academic) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = target)) +
  #geom_point(aes(group = target)) + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  coord_cartesian(ylim = c(-.15, .3)) +
  scale_color_manual(values = pal_self_other) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  coord_cartesian(ylim = c(-.15, .3)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  #geom_point() + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))

self/other evaluation task analyses

visualize raw

age

(status_plot_raw_age = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(age, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(status_plot_raw = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(pdss, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(task_percent ~ subjectID, data = status_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)
## [1] 0.18

age

`status model 1` = lmer(task_percent ~ 1 + target*domain + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 2 age` = lmer(task_percent ~ 1 + target*domain + age_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 age` = lmer(task_percent ~ 1 + target*domain + age_c + age_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            target*domain*age_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`status model 2 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + pdss_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + target*domain*pdss_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`status model 1`, `status model 2 age`, `status model 3 age`, `status model 4 age`, `status model 5 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
status model 1 15 5237.97 – – –
status model 2 age 16 5239.91 0.05 1 .815
status model 3 age 17 5241.69 0.22 1 .637
status model 4 age 19 5233.46 12.23 2 .002
status model 5 age 23 5231.07 10.39 4 .034

puberty

anova(`status model 1`, `status model 2 puberty`, `status model 3 puberty`, `status model 4 puberty`, `status model 5 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
status model 1 15 5237.97 – – –
status model 2 puberty 16 5238.96 1 1 .317
status model 3 puberty 17 5238.77 2.2 1 .138
status model 4 puberty 19 5232.81 9.95 2 .007
status model 5 puberty 23 5237.02 3.8 4 .434

summarize best fitting models

age

`status model 5 age` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, age 13)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (academic, other, age 13) 57.87 [51.28, 64.47] 3.36 17.21 174.40 < .001
Target (self) 4.82 [-4.34, 13.97] 4.67 1.03 182.28 .304
Domain (social) 2.29 [-6.42, 11.01] 4.45 0.52 188.62 .607
Age -0.56 [-1.65, 0.53] 0.55 -1.01 365.43 .312
Age2 0.14 [-0.42, 0.69] 0.28 0.48 347.30 .631
Target (self) x Domain (social) 18.51 [6.41, 30.61] 6.17 3.00 194.63 .003
Target (self) x Age -0.68 [-2.22, 0.86] 0.79 -0.87 366.83 .386
Domain (social) x Age 0.55 [-0.97, 2.07] 0.78 0.70 371.77 .481
Target (self) x Age2 0.43 [-0.35, 1.21] 0.40 1.07 348.69 .283
Domain (social) x Age2 -0.48 [-1.26, 0.29] 0.40 -1.22 352.91 .222
Target (self) x Domain (social) x Age 2.09 [-0.04, 4.22] 1.09 1.92 372.86 .056
Target (self) x Domain (social) x Age2 -0.64 [-1.73, 0.46] 0.56 -1.14 356.85 .255
print(VarCorr(`status model 5 age`), comp = c("Variance","Std.Dev."), digits = 5)
##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              377.05  19.418                       
##            targetself               689.50  26.258   -0.815              
##            domainsocial             565.08  23.771   -0.866  0.798       
##            targetself:domainsocial 1040.78  32.261    0.772 -0.936 -0.905
##  Residual                           239.13  15.464

puberty

`status model 4 puberty` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, stage 3)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", term),
         term = gsub("pdss_c", "Puberty", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (academic, other, stage 3) 59.23 [54.19, 64.28] 2.57 23.01 78.80 < .001
Target (self) 8.48 [1.51, 15.45] 3.56 2.38 80.81 .02
Domain (social) -1.54 [-8.03, 4.94] 3.31 -0.47 80.29 .642
Puberty -1.68 [-3.91, 0.54] 1.13 -1.48 380.92 .139
Target (self) x Domain (social) 12.33 [3.43, 21.24] 4.54 2.71 80.80 .008
Target (self) x Puberty -1.14 [-4.28, 2.00] 1.60 -0.71 379.97 .478
Domain (social) x Puberty 1.37 [-1.72, 4.46] 1.57 0.87 382.90 .385
Target (self) x Domain (social) x Puberty 3.53 [-0.78, 7.83] 2.20 1.60 383.96 .109
print(VarCorr(`status model 4 puberty`), comp = c("Variance","Std.Dev."), digits = 5)
##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              367.36  19.167                       
##            targetself               687.68  26.224   -0.825              
##            domainsocial             561.36  23.693   -0.869  0.796       
##            targetself:domainsocial 1030.63  32.103    0.797 -0.941 -0.908
##  Residual                           244.39  15.633

visualize fitted

age

model_to_plot = `status model 5 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_age = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(age, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = "none"))

puberty

model_to_plot = `status model 4 puberty`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_pdss = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(pdss, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = c(.8, .14),
          legend.spacing.x = unit(.1, 'cm'),
          legend.box = "horizontal"))

combined

cowplot::plot_grid(status_plot_age, status_plot_pdss, labels = c('A', 'B'), ncol = 2)

SPPC competence analysis

visualize raw

age

(comp_plot_raw_age = comp_model_data %>%
  ggplot(aes(age, competence, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(comp_plot_raw = comp_model_data %>%
  ggplot(aes(pdss, competence, color = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(competence ~ subjectID, data = comp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)
## [1] 0.31

age

`comp model 1` = lmer(competence ~ 1 + domain + 
                              (1 | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 2 age` = lmer(competence ~ 1 + domain + age_c + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 age` = lmer(competence ~ 1 + domain + age_c + age_c2 + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 age` = lmer(competence ~ 1 + domain*age_c + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 age` = lmer(competence ~ 1 + domain*age_c + domain*age_c2 + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`comp model 2 puberty` = lmer(competence ~ 1 + domain + pdss_c + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 puberty` = lmer(competence ~ 1 + domain + pdss_c + pdss_c2 + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 puberty` = lmer(competence ~ 1 + domain*pdss_c + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 puberty` = lmer(competence ~ 1 + domain*pdss_c + domain*pdss_c2 + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`comp model 1`, `comp model 2 age`, `comp model 3 age`, `comp model 4 age`, `comp model 5 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model model df AIC X2 X2 df p
comp model 1 4 430.67 – – –
comp model 2 age 7 429.81 6.85 3 .077
comp model 3 age 8 428.13 3.68 1 .055
comp model 4 age 8 431.13 0 0 –
comp model 5 age 10 431.05 4.08 2 .13

puberty

anova(`comp model 1`, `comp model 2 puberty`, `comp model 3 puberty`, `comp model 4 puberty`, `comp model 5 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model model df AIC X2 X2 df p
comp model 1 4 430.67 – – –
comp model 2 puberty 7 431.15 5.51 3 .138
comp model 3 puberty 8 432.13 1.03 1 .311
comp model 4 puberty 8 432.73 0 0 –
comp model 5 puberty 10 435.07 1.65 2 .437

summarize best fitting model

age

# summarize best fitting model
`comp model 3 age` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (social)", term),
         term = gsub("age_c", "age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--"),
         SE = sprintf("%.3f", SE)) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
intercept (academic, age 13) 3.25 [3.16, 3.34] 0.046 71.20 249.20 < .001
domain (social) -0.22 [-0.30, -0.14] 0.041 -5.38 228.29 < .001
age 0.02 [-0.00, 0.03] 0.010 1.63 73.36 .106
age2 -0.01 [-0.02, 0.00] 0.004 -1.92 315.15 .056
# summarize random effects
print(VarCorr(`comp model 3 age`), comp = c("Variance","Std.Dev."), digits = 4)
##  Groups    Name        Variance Std.Dev. Corr 
##  subjectID (Intercept) 0.022106 0.14868       
##            age_c       0.001837 0.04286  -0.21
##  Residual              0.152041 0.38992

visualize fitted

age

model_to_plot = `comp model 3 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
comp_model_data$expected = predict(model_to_plot, newdata = comp_model_data, re.form = reForm)
comp_model_data$expected_mean = predict(model_to_plot, newdata = comp_model_data, re.form = NA)

(comp_plot_age = comp_model_data %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "SPPC competence rating\n") +
    dcbw +
    theme(legend.position = "none"))

SPPC importance analysis

visualize raw

age

(imp_plot_raw_age = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(age, importance, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    labs(x = "\npubertal stage (PSD score)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(imp_plot_raw = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(pdss, importance, color = domain)) +
    geom_jitter(alpha = .1, width = .05, height = .05) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(importance ~ subjectID, data = imp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)
## [1] 0.24

age

`imp model 1` = lmer(importance ~ 1 + domain +
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 2 age` = lmer(importance ~ 1 + domain + age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 age` = lmer(importance ~ 1 + domain*age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`imp model 2 puberty` = lmer(importance ~ 1 + domain + pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 puberty` = lmer(importance ~ 1 + domain*pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`imp model 1`, `imp model 2 age`, `imp model 3 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
imp model 1 6 356.28 – – –
imp model 2 age 7 356.65 1.62 1 .203
imp model 3 age 8 356.69 1.97 1 .161

puberty

anova(`imp model 1`, `imp model 2 puberty`, `imp model 3 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
imp model 1 6 356.28 – – –
imp model 2 puberty 7 357.77 0.51 1 .476
imp model 3 puberty 8 359.27 0.5 1 .482

summarize best fiting model

`imp model 1` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (social)", term),
         term = gsub("age_c", "age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
intercept (academic, age 13) 3.35 [3.23, 3.48] 0.06 52.25 55.06 < .001
domain (social) -1.01 [-1.21, -0.81] 0.10 -10.00 55.63 < .001
# summarize random effects
print(VarCorr(`imp model 1`), comp = c("Variance","Std.Dev."), digits = 4)
##  Groups    Name         Variance Std.Dev. Corr 
##  subjectID (Intercept)  0.04644  0.2155        
##            domainsocial 0.19984  0.4470   -0.24
##  Residual               0.29299  0.5413

visualize fitted

imp_model_data_na = imp_model_data %>%
  filter(!is.na(importance))

model_to_plot = `imp model 1`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
imp_model_data_na$expected_mean = predict(model_to_plot, newdata = filter(imp_model_data_na, !is.na(importance)), re.form = NA)

(imp_plot_age = imp_model_data_na %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "SPPC importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .14),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

SPPC competence importance plot

cowplot::plot_grid(comp_plot_age, imp_plot_age, ncol = 2, labels = c("A", "B"))